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Abstract. We present sCOLA - an extension of the N-body COmoving Lagrangian Accelera¬ 
tion (COLA) method to the spatial domain. Similar to the original temporal-domain COLA, 
sCOLA is an N-body method for solving for large-scale structure in a frame that is comoving 
with observers following trajectories calculated in Lagrangian Perturbation Theory. Incor¬ 
porating the sCOLA method in an N-body code allows one to gain computational speed by 
capturing the gravitational potential from the far field using perturbative techniques, while 
letting the N-body code solve only for the near field. The far and near fields are completely 
decoupled, effectively localizing gravity for the N-body side of the code. Thus, running an 
N-body code for a small simulation volume using sCOLA can reproduce the results of a stan¬ 
dard N-body run for the same small volume embedded inside a much larger simulation. We 
demonstrate that sCOLA can be safely combined with the original temporal-domain COLA. 
sCOLA can be used as a method for performing zoom-in simulations. It also allows N-body 
codes to be made embarrassingly parallel, thus allowing for efficiently tiling a volume of in¬ 
terest using grid computing. Moreover, sCOLA can be useful for cheaply generating large 
ensembles of accurate mock halo catalogs required to study galaxy clustering. Surveys that 
will benefit the most are ones with large aspect ratios, such as pencil-beam surveys, where 
sCOLA can easily capture the effects of large-scale transverse modes without the need to 
substantially increase the simulated volume. As an illustration of the method, we present 
proof-of-concept zoom-in simulations using a freely available sCOLA-based N-body code. 


1 Introduction 


Understanding Large Scale Structure (LSS) of the universe has been a multi-generational ef¬ 
fort in understanding the macroscopic physics governing the evolution of the universe, arising 
from a trivial microscopic law - the law of universal gravitation (corrected for the expansion 
of the universe, as well as baryonic physics at small scales). Analytical approaches to the 
problem have recently (e.g. [1-3]) produced new insights into the physics of LSS, signifi¬ 
cantly improving the predictive power of analytical techniques. Yet, numerical simulations 
are indispensable when one needs to go beyond the mildly non-linear regime, as self-consistent 
analytical approaches break down beyond the non-linear scale, corresponding to a wavevector 
of ^ 0.5Mpc jh at redshift of zero. 

As the volume and precision of observational surveys of LSS increase, modelling LSS 
observables (e.g. galaxies, weak lensing) has been pushing the limits of both analytical and 
numerical techniques. One of the big challenges is accurately modelling uncertainties. That 
requires reducing sample variance for the 4-point function (covariance matrices) of the tracer 
field of LSS of choice for a particular survey. The problem can in principle be solved by 
running hundreds and even thousands of numerical simulations, but until recently this was 
feasible only after introducing rather crude approximations (e.g. the PTHalos method [4, 5]). 

The push to construct crude, yet realistic, numerical simulations to tackle the problem 
of modeling covariances has led to the development of numerous prescriptions [6-12] for fast 
generation of density fields and/or mock halo and galaxy catalogs (for a recent review, see 
[13]). One of them, the COmoving Lagrangian Acceleration (COLA) method ([6], hereafter 
referred to as TZE), is a proposed modification to existing N-body codes, which allows for 
the reuse of post-processing pipelines, designed to work with N-body codes. COLA allows for 
a fine control over the trade-off between accuracy and speed and as a cheap N-body method 
it has been successfully applied in extracting cosmological information from observations 
[14-18]. 

In its original formulation, the COLA method used the fact that the temporal evolution 
of large-scale modes is well described by Lagrangian Perturbation Theory (LPT) [19] (see 
[20] for a recent detailed discussion of LPT). COLA decouples the temporal evolution of the 
large and small scales in N-body codes by evolving the large scales using the exact LPT 
results for the growth factor of those scales, while letting the N-body code calculate the time 
evolution of the small scales. 

In its approach of combining Perturbation Theory (PT) with N-body simulations, COLA 
is taking advantage only of the fact that the temporal behavior of large scales can be modeled 
by LPT. Yet, there is still room for improvement, as PT can also be used to describe how 
distant objects gravitatonally influence local dynamics. Such far-held effects are dominated 
by large scales - most importantly by rescaling the local matter density and introducing tidal 
effects. The reason for that is that the detailed small-scale distribution of matter far away 
from a region of interest only changes the higher multipoles (in a multipole expansion) of the 
potential which decay much faster than the low multipoles (which are affected by the large 
scales). 

A method that allows for the complete decoupling of local dynamics (modeled with 
an N-body) from far-held effects (modeled with PT) can have numerous applications. An 
obvious one is performing zoom-in simulations without the need to substantially extend the 
simulation volume around a region of interest. It would also allow N-body codes to be made 
embarrassingly parallel, thus allowing for efficiently tiling a volume of interest. Moreover, 
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such a method can be useful for cheaply generating large ensembles of accurate mock halo 
catalogs required to study galaxy clustering. Surveys that would benefit the most are ones 
with large aspect ratios, such as pencil-beam surveys, where such a method can capture the 
effects of large-scale transverse modes in PT without the need to substantially increase the 
simulated volume 1 . 

In this proof-of-concept paper we propose exactly such a modification to existing N- 
body codes by extending the COLA method to the spatial domain, resulting in what we call 
spatial COLA (sCOLA). In Section 2 we present the general idea behind the sCOLA method. 
We present results from a set of zoom-in sCOLA-based N-body simulations in Section 3, and 
we summarize in Section 4. We leave the detailed description of the sCOLA method to the 
appendices. 

2 Augmenting Simulations with Perturbation Theory to Model the Far 
Field 

2.1 Review of the original temporal COLA method 

Since sCOLA is an extension of COLA, we start with a quick review of the original COLA 
method. In COLA (see TZE) one starts with the standard equation of motion for cold dark 
matter (CDM), which schematically states 

d 2 t x = -VV“ 2 5[*](*) . (2.1) 

In COLA, one rewrites the above equation as 

c? 2 x res = —VV _2 5[x](x) — 5 2 xlpt , with x res = x — xlpt • (2.2) 

In the above, VV _2 5[a](&) is a short-hand for the gradient of the inverse Laplacian of the 
fractional CDM overdensity arising from a distribution of particles at positions a, evaluated at 
position b of the particle for which the equation of motion is solved for. Here x(g, t) = xlpt + 
x res is the Eulerian position of CDM particles at time £, which started out at (Lagrangian) 
position q\ xlpt is the LPT approximation to x; and x res is the (residual) displacement of 
CDM particles as measured in a frame comoving with “LPT observers”, whose trajectories 
are given by xlpt- Note that xlpt(<Z,£) is the perturbative solution to (2.1): 

d 2 xppT = - VV 2 ^[^lpt](^lpt) (solved in PT), (2-3) 

where the equation above is solved using the standard LPT prescription of expanding it in 
powers of the linear overdensity. 

In COLA, one discretizes the d 2 operator on the left-hand side of (2.2), while calculating 
dt x lpt on the right-hand side exactly in LPT. Therefore, for zero timesteps, COLA recovers 
the LPT results, which implies that large scales are followed with an accuracy of at least 
that of LPT. Corrections to LPT are obtained by increasing the number of timesteps - a 
characteristic especially important for recovering the correct time evolution of small scales. 
When the number of timesteps becomes large, COLA reduces to a standard N-body code. In 
the intermediate regime, (for 0(10) timesteps) COLA is able to give a good approximation 
of halo statistics at the expense of not resolving the detailed particle trajectories inside halos 
(see TZE and [13] for further discussion). 

An approach to solve this problem has been put forward by [21], which uses the properties of the periodic 
boundary conditions to map a simulation done on a periodic cube into elongated box-like shapes. However, this 
approach breaks the behavior of the large modes by changing the topological connectedness of the simulated 
region. 
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2.2 Introducing sCOLA 

As discussed in the introduction, the sCOLA method is a modification to existing N-body 
codes allowing one to perform simulations without the need to substantially extend the 
simulated volume beyond the region of interest in order to capture far-held effects, such 
as the effects of super-box modes (modes longer than the size of the simulation box). To 
understand the principles behind sCOLA, we will first state the equations behind sCOLA 
and will only then proceed to discuss their interpretation and characteristics in detail. 

Let us imagine that we want to simulate the detailed dynamics of a region of interest 
(referred to as the COLA volume/box from now on), which is embedded in a much larger 
region (the “full” volume/box). sCOLA stipulates that one can then model the dynamics 
inside the COLA volume by rewriting the equation of motion (2.1) as (we reintroduce the 
omitted constants and Hubble expansion in Appendix A): 

Xres — ^COLA^K*) - ^ L C p°t LA , with ^res — #lpt ^ (2.4) 

where we introduced an auxiliary LPT displacement field, #lpt jA (< 7? which similar to (2.3) 
is defined as the perturbative solution to (we reintroduce the omitted constants and Hubble 
expansion in Appendix B) 

^£p°t LA = -VVc 2 OL a^lpt](*lpt) (solved in PT). (2.5) 

The subscript COLA in Vqqla stands for the fact that the inverse Laplacian used to obtain 
the gravitational potential from the density field is restricted only over the near field (the 
COLA volume). In other words, the potential computed with Vqq LA 5 includes only particles 
in the COLA volume, not the full volume. Compare that to V -2 in (2.2), which is over the 
full box. This marks the crucial difference between the original temporal domain COLA and 
sCOLA. Equations (2.4) and (2.5) solved for x for the particles inside the COLA volume 
summarize the sCOLA method. 

One can think of the sCOLA equations in the following way. We calculate the acceler¬ 
ation of each particle in the COLA box in the frame of an accelerated “LPT observer”, who 
follows a trajectory (^lpt^) as specified by LPT in that same COLA volume. We then use 
that acceleration to calculate the extra displacement (x res ) of the particle, which we need to 
add to the LPT displacement in order to find the actual position of the particle. However, 
instead of adding x res to ^lpt^’ which is calculated in the COLA box, we add x res to the 
LPT position (xlpt) as calculated in the full box. 

Simply put, the idea is that by using the wrong observer (^lpt" A ) i n the wrong universe 
(by calculating accelerations in the COLA box) and then pretending we did no such thing 
(by adding the true LPT position, cclpt), we make two mistakes, which will tend to cancel 
out. We demonstrate that that is indeed what happens in the rest of this paper. 

2.3 Discussion 

There are several important aspects of the sCOLA equations that we highlight below 2 : 

1) Once the initial conditions xlpt are known, the sCOLA equations above make no 
reference to particles or fields (be that density, potential or otherwise) outside of the COLA 

2 Note that sCOLA is somewhat similar to existing post-processing methods for capturing the effects on 
N-body simulations from changes in the background cosmology [22] or for resampling the large-scale modes of 
existing simulations [23]. However, unlike sCOLA, both of those methods fail to capture the coupling between 
the large- and short-scale displacement fields and we will, therefore, not consider them further. 
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volume. Thus, the non-locality of gravity arising from the inverse Laplacian in the calculation 
of the gravitational potential, has been severely limited to the COLA volume. Therefore, 
a sCOLA-based N-body code solving for the positions of particles inside the COLA volume 
need only work with particles and fields inside of that small volume and not the full volume, 
which the COLA volume is embedded in. Thus, the main objective of sCOLA - to decouple 
the near and far fields, using PT for the far field, and an N-body for the near field - has been 
achieved. 

2) As mentioned above, the only place where the far field enters is when one calculates 
the initial conditions xlpt i n the full box, as the LPT displacements of the particles in the 
COLA box are correlated with the LPT displacements of all particles in the full box (albeit 
to different degrees). One should generate those initial conditions in the same way one would 
for a standard N-body run. For zoom-in simulation, for example, one can use multi-scale 
initial conditions generators such as those presented in [24] or [25]. 

3) The sCOLA equations become exact in the limit when the COLA volume matches 
the full volume, reducing to the original COLA equation, eq. (2.2). When those two volumes 
do not coincide, the error one makes is given by the residual acceleration defined by: 


ares(q) = ^VV cola 5[xlpt](^LPt) - VV 5[xlpt](#LPt)J 

- (w^ola^WN-VV-^M^)) , 


in PT 


(2.6) 


where the first term is calculated in PT, and the right-hand side is evaluated for the parti¬ 
cle/LPT observer starting at q. Clearly, if the COLA volume matches the full volume, then 
the inverse Laplacians coincide, a res vanishes, and sCOLA reduces to the original COLA 
method as mentioned above. When the COLA volume does not cover the full volume, then 
a res represents corrections coming only from the far field (the field not covered by the COLA 
volume). As LPT is an extremely good approximation to the density field at large scales 
(e.g. [20]), then a res inside the COLA volume contains only corrections from the far field that 
come only from small-scale deviations of the density field from its LPT prediction. Thus, as 
long as the COLA box is not too small (see below), then a res will contain only the gradients 
of quickly decaying high-order multipoles of the far field. 

4) From the discussion in (3) above, we see that the errors arising from a res are most 
prominent near the boundaries of the COLA volume. Thus, one should always make sure 
to pad their region of interest by roughly the scale below which LPT breaks down (see also 
(3) above), which is about lOMpc jh at redshift zero (e.g. [26]). One should also make sure 
that their COLA volume encloses all particles that eventually make it into their region of 
interest. The typical distance CDM particles travel to redshift zero is again about lOMpc jh 
(e.g. [20]). Thus, as a rule of thumb, for a simulation down to redshift of zero, the COLA 
volume should cover a region which is roughly ^lOMpc /h larger in all directions around the 
region of interest. The actual number will depend on the specific application of sCOLA and 
the desired accuracy. 

5) One may wonder what boundary conditions (BCs) one should use for the COLA 
volume. The short answer is: any will probably do as long as those BCs are applied con¬ 
sistently in both sCOLA equations. For any reasonable COLA volume (see above), errors 
arising from the BCs will only modify the amplitudes of the far-field high-order multipoles 
entering a res , which, as we already discussed, decay rapidly away from the boundaries of the 
COLA volume. 
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Lcola 

[Mpc//i] 

nt 

Mass 
[M q /h] 

Mean 

Position 
[Mpc /h\ 

Size 

[Mpc /h\ 

Mean 
Velocity 
[Mpc /h\ 

Velocity 
Dispersion 
[Mpc /h\ 

Reference results (L =Lcola, nt — 100) 

100 

100 

2.08xl0 14 

0 

0.79 3.18 

9.65 

Original temporal COLA results relative to reference 

L =Lcola 

100 

15 

-15 % 

0.10 

+15 % 

0.23 

-17 % 

100 

10 

-18 % 

0.17 

+25 % 

0.30 

-21 % 

100 

7 

-38 % 

0.06 

+31 % 

0.43 

-30 % 

sCOLA results relative to reference 
n t = 100 

75 

100 

-1 % 

0.01 

+1 % 

0.06 

-0.3 % 

50 

100 

-0.1 % 

0.03 

+2 % 

0.04 

-0.4 % 

25 

100 

-2 % 

0.25 

+68 % 

0.37 

-11 % 

sCOLA results relative to reference 
nt = 10 

100 

10 

-18 % 

0.17 

+25 % 

0.30 

-21 % 

75 

10 

-19 % 

0.18 

+26 % 

0.30 

-22 % 

50 

10 

-19 % 

0.18 

+26 % 

0.25 

-22 % 

25 

10 

-20 % 

0.50 

+30 % 

0.20 

-24 % 


Table 1: The table summarizes the properties of the central halo as obtained from each 
snapshot (identified by nt and Lqola) shown in the figures. The halo properties shown are 
the halo mass, the center-of-mass position, the halo size (defined as the root-mean-square 
halo particle positions); the center-of-mass halo velocity, scaled to show the shift of the 
halo position in redshift space; and the velocity dispersion (defined as the rms halo particle 
velocity). Apart from the reference results, all other results for the halo properties show the 
error with respect to the reference result. The errors are absolute unless shown as percentage 
difference. Reducing the number of timesteps for COLA naturally degrades the quality of 
the simulation. However, note that the sCOLA-side of the method proves robust to reducing 
the size of the COLA box until the COLA boundary region (at redshift zero, corresponding 
to the region ^ lOMpc jh or less from the edge) invades the region of interest. 


6) The time operator in the sCOLA equations can be discretized either in the standard 
approach (discretize all <9|), or in the original temporal domain COLA approach, where one 
discretizes only on the left-hand side of (2.4) as the time evolution of xlpt and ^lpt^ 
known to be governed by the usual linear (and second order, for second order LPT) growth 
factor. In this paper we use the latter approach which will allow us to demonstrate that 
the temporal and spatial versions of COLA can be used successfully both separately and in 
conjunction. 
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Figure 1: Slices of thickness 12.5Mpc jh at z — 0 through one and the same CDM density 
field as obtained with the original time-domain-only COLA method by varying the number 
of timesteps, n*, across columns. Rows show the result at different magnifications, with 
the top row slices being 25Mpc jh on the side, and the bottom row slices being 12.5Mpc jh 
on the side. One should note that for rtt < 10 halos puff up significantly, rendering their 
identification more unreliable. We use the same color map for all projected density snapshots 
in this paper. 


3 Illustrative example 

To illustrate how sCOLA can be implemented in an N-body code, we developed a freely 
available and open-source parallel Python 3 code, implementing in conjunction the sCOLA 
and the original COLA methods, dubbed pyCOLA. To generate the initial conditions for 
the simulations presented below, we use the multi-scale initial conditions generator, MUSIC 4 
[24]. pyCOLA then calculates the LPT displacements in the COLA volume, #lpt" A ’ by 
implementing the equations in Appendix B, corresponding to the schematic eq. (2.5). py¬ 
COLA then solves for the particle trajectories in the COLA box by evolving the particle 
positions and velocities using the kick-drift-kick method in Appendix A.3, corresponding to 
the schematic eq. (2.4). pyCOLA implements periodic BCs for the COLA volume and the 
force calculation is done using a PM grid over the same COLA volume. We only include 
particles that are present in the COLA volume in the initial conditions, and do not include 
particles that eventually enter that region from the outside. 

The second order initial conditions (giving xlpt) in the full box are calculated 5 with 
MUSIC using 512 3 particles in the full box of size L = 100Mpc//i. The central 45% 

3 https://www.python.org/ 

4 Note that while working on this paper, we encountered a bug in MUSIC, which resulted in a wrong 
second-order displacement field. O. Hahn has kindly found and fixed the bug in version 1.51 of MUSIC. 

5 The adopted cosmology in standard notation is: Pm = 0.274, Pa = 0.726, P b — 0.046, h = 0.7, as = 0.8, 
n s = 0.95. 


- 7 - 













I-.V.A'- lou \I:»> ; v ’■' > .4 4 L C ola= 75 Mpc/?» 


1 m; ■■ t 


> . 


%>Q 




//•,..•-->100 

CopA = 100Mpc//i ; 




72a =100 


LcopA ^T5. -Mpc/' h ; 




'$t « 


72^=100 


.'. jr 


- *;* • —«j- ! -A——rr"—- . — 

^coi/A Mpe//i 4|j ^coi/a = AP Mp &lh 


^^ 74 ... 


■■$&$■' . V;. 




= 100 


^COLA - 50 Mpc/fl 

■fr ^ 

r-. 


*£*£# L 


72a =100 


Lcola — 50 Mpc/fi 


;• v’\. . v -; 

| i 

i ;■ 


n, =100 


^COLA - 25 Mpc/fi 


72a =100 


Lcola — 25 Mpc/fi 



72a =100 


i'ropA — ip9 Mpc fh W. ;• 


^cola - 25 Mpc/fi 





..• yV-\ • %; rV-.y y . . §■?• \ /•’, . * 

v- :;: v r $&> -V ;>y i >- 


--'.V'" I 

■ K 

:#y^.roo 


a ■ a : : 

. .. i*r^-ioo -*■, ■•' 

Jj , ?*. : . -. • .,. »j ,, 

:pc/ftV A ' . i.,,„ * Tr, Nrprfc V. A * v > 

'*V. \ •' • , - ’ ( - r "*V^' : " 



Lcola= 100 Mpem 


V 




l coi'a= 50 M P.c//i v ,1 * > L cola =25 Mpc/ h 



Figure 2: Slices of thickness 12.5Mpc jh at z = 0 through the same CDM density field 
as obtained with the proposed extension of the COLA method to the spatial domain. The 
number of timesteps for all slices is fixed at nt = 100. Columns show different choices of 
Lcola- Rows show the result at consecutive magnifications, with the slices being 100 (top 
row), 75, 50, 25, 12.5 (bottom row) Mpc jh on the side. For Lcola = 25Mpc//i, not all 
particles in the vicinity of the central halo are present in the initial Lagrangian volume, 
which leads to noticeable differences in the distribution of substructures. Note that the slices 
in the last two rows are offset relative to the other rows, so that they are centered around 
the most massive halo. 


































(lengthwise) of the full box are refined with a grid of particles with a particle spacing of 
L/1024 = 0.098 Mpc//i, corresponding to particle masses of 7.08 x 10 7 M Q //i. We evolved 
those same initial conditions with pyCOLA from redshift of 9 down to redshift of zero us¬ 
ing various number of timesteps, n^, and COLA box sizes, Lcola- The LPT displace¬ 
ment field inside the COLA volume was calculated to second order with the 2-pt rules in 
Appendix B using pyCOLA. We then evolved the particle positions using pyCOLA with 
nt — 100,15,10, 7 and Lcola = 100, 75, 50, 25 Mpc jh. The simulations were centered around 
a halo of mass 2.08 x 1O 14 M 0 //i. We used a fixed force resolution (PM grid spacing) of 
L/1536 = 0.065 Mpc//i for all runs. Note that both the size of the PM grid and the number 
of particles inside the COLA volume change with Lcola as we are keeping the force reso¬ 
lution and particle masses constant, while changing the volume (corresponding to Lcola) 
simulated by the N-body code. 

In Figure 1 we show slices through the CDM density field obtained with the original 
temporal COLA method 6 by setting (see Section 2.3) the full box and COLA box sizes to be 
the same, i.e. L =Lcola = 100 Mpc /h. We vary the number of timesteps between 100 and 
7 across the columns of that panel plot. The two rows correspond to two different magnifi¬ 
cations: the top slices are 25 Mpc jh on the side, while the bottom are 12.5 Mpc jh. One can 
clearly see that the original temporal COLA method behaves as expected. It reproduces the 
large-scale features faithfully, while performing gradually worse at small scales when n t is 
lowered. For nt < 10 halos puff up significantly as the few timesteps become insufficient for 
the N-body side of the code to keep particles tightly bound to halos (see TZE). One can see 
that from Table 1 (discussed in more detail below) as well, where the central halo mass shows 
a gradual decline with decreasing nt, while the halo size increases. Note that the errors in 
the halo position in both real and redshift space are within 500 kpc//i, while the finger-of-god 
(FoG) effect is consistently underestimated. For a discussion of how these imperfections can 
be fixed, see TZE for example. 

Next we focus on the simulation runs using the new sCOLA method. The results 
for nt = 100 are presented in Figure 2, while those for nt = 10 - in Figure 3. Columns 
represent different choices of COLA volume size, Lcola > while rows represent consecutive 
magnifications. Concentrating on the first row of each figure, one can see how the four 
different simulation volumes compare in size. Inspecting the last row, one can see that the 
halo substructure is nearly identical across all panels with the exception of the rightmost panel 
(see below). Concentrating on the intermediate rows, we can see that large-scale structure is 
captured faithfully by sCOLA, and no large-scale distortions are introduced in the bulk even 
using clearly suboptimal periodic boundary conditions for the COLA volumes. 

The rightmost column (Lcola = 25 Mpc /h) shows clear problems: some substructure 
is missing, and the density field near the boundaries of the COLA volume is too smooth. 
The reason is straightforward to pinpoint and was already anticipated in Section 2.3. The 
COLA box for that run is too small to include all particles that end up in the vicinity of the 
central halo, thus providing an illustration of how sCOLA behaves when it is pushed beyond 
breaking (see Section 2.3). 

It is not surprising that the snapshot with Lcola = 25 Mpc /h exhibits problems. The 
boundary region, which is the region within ^ 10 Mpc jh from the COLA boundary, where 
we expect to see sCOLA misbehaving (see Section 2.3), occupies most of the volume in this 


6 Note that we did not try to optimize the choice for tilpt for this simulation set-up as prescribed in TZE, 
so the results we present here using the original COLA method can be suboptimal. 
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Figure 3: The same as in Fig. 2 but for = 10. This figure illustrates that COLA can be 
safely applied in the spatial and temporal domains simultaneously. 
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snapshot. Thus, it is only for Lcola ^ 50 Mpc /h and above that we can expect the boundary 
effects on the central halo to be small, which is indeed confirmed. 

Another example of sCOLA breaking is close to the COLA volume boundaries, which 
we again anticipated in Section 2.3. In Figures 2 and 3, one can clearly see artifacts from 
the non-zero a res and the periodic boundary conditions. Those are especially prominent on 
the right edge of the Lcola = 50 Mpc jh COLA volume, where one can see a high-density 
wall of particles, which is missing from the larger COLA volume snapshots. Inspecting the 
larger COLA volume snapshots, we see that those particles are actually flowing to the left, 
outside of the central 50 Mpc//i, and are therefore wrapped around by the periodic BCs in the 
Lcola = 50 Mpc jh snapshot. Similarly, the lower COLA boundary for the same snapshot 
is clearly emptied out since in our pyCOLA implementation no particles from outside the 
COLA volume are allowed to flow in. 

Nevertheless, those boundary artifacts in the Lcola = 50 Mpc jh snapshot do not affect 
the accurate recovery of the density field in the bulk. Thus, in this example sCOLA has 
allowed us to reduce the simulation volume by a factor of 8 without sacrifice in the recovered 
central halo bulk properties, environment or substructures. 

For nt — 100 the temporal domain behavior of sCOLA is approaching a full-blown 
N-body simulation. We next show how sCOLA behaves when nt is substantially reduced, 
thus making full use of the original temporal-domain COLA side of the method. Figure 3 
demonstrates this feature of the method by significantly reducing both Lcola and nt. As 
one can see, those combined sCOLA/COLA snapshots share the same pros and cons of the 
two methods. Thus, depending on the application, one can apply COLA independently in 
the spatial and temporal domains, or simultaneously, as needed. 

To be more quantitative, we investigate the properties of the central halo. We ran a 
friends-of-friends halo finding algorithm 7 on all simulation snapshots presented in this paper, 
using a fixed linking length of 0.2. The results for the halo at the center of each snapshot 
are collected in Table 1. As reference values, we use the results from the most accurate 
simulation, corresponding to L =Lcola = 100 Mpc //i, nt = 100. Thus, the reference value 
for the mass of the central halo is 2.08 x 10 14 Mq//i. In the table we list the error in the 
recovered halo mass relative to the reference value, the error in its center-of-mass position and 
velocity 8 * * * and the error in the root-mean-square (rms) of the halo particle positions, giving 
a measure of the halo size, as well as the error in the rms halo particle velocities, giving a 
measure of the FoG effect in redshift space. 

In Table 1, we can see that for nt = 100, sCOLA performs superbly for Lcola = 
75 Mpc jh and Lcola = 50 Mpc//i, i.e. as long as all halo particles are captured inside the 
COLA volume. All errors in both real and redshift space are around or well within the force 
resolution, and the halo mass is recovered to 1% or better. 

When nt is pushed down to 10, we see that most of the deterioration can be attributed to 
the smaller number of timesteps (i.e. to the original COLA method), and not to the sCOLA 
side of the code. That is evident when comparing the results of smaller COLA volumes to 
the result of the Lcola = 100 Mpc /h, nt = 10 snapshot. 

7 We use the University of Washington FOF code available here: http://www-hpcc. astro .Washington, 
edu/tools/fof.html 

8 For velocity, we quote values for (dx / drf) / (aH (a)), where H(a) is the Hubble parameter at scale factor 

a, 7] is conformal time, and x is comoving distance. Thus, the units of our velocity variable are Mpc /h. This 

definition allows calculating redshift-space positions trivially: one simply has to add the line-of-sight velocity 

to the halo position. 
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4 Summary 


In this proof-of-concept paper we presented sCOLA - an extension to the N-body COLA 
method in the spatial domain. sCOLA is a straightforward modification to N-body codes, 
which uses perturbation theory to calculate the far field potential from masses outside a region 
of interest, while solving for the near field using an N-body code. The use of perturbation 
theory allows for the N-body side of sCOLA to perform all calculations entirely within the 
region of interest (dubbed the “COLA volume”), without the need to substantially increase 
the simulated volume to capture the effects of large modes or structures outside that region. 

Compared to standard zoom-in simulation methods, the main advantage of sCOLA is 
that it makes the different volumes completely independent once the initial conditions are 
computed. Thus, there is no need to add lower resolution zones around the region of interest 
in sCOLA. 

The sCOLA method leads to the exact equations of motion when the COLA volume 
coincides with the full volume, capturing all large-scale modes affecting the region of interest. 
When the COLA volume does not cover the full volume, the approximations behind sCOLA 
require one to run the N-body code behind sCOLA on a region which is roughly lOMpc jh 
larger in all directions than the region of interest. 

We demonstrated the usage of sCOLA by showing results from zoom-in simulations 
performed with a freely available sCOLA-based N-body code (called pyCOLA), which we 
developed for this paper. Our results demonstrate that sCOLA faithfully recovers large-scale 
structure inside a COLA volume as small as 50Mpc jh (possibly even smaller), while keeping 
substructures intact. Moreover, we showed that sCOLA can be safely used in conjunction 
with the original COLA method. We find that reducing the number of timesteps for COLA 
naturally degrades the quality of the simulation. However, the sCOLA-side of the method 
proves robust to reducing the size of the COLA box until the COLA boundary region (at 
redshift zero, corresponding to the region ~ lOMpc /h or less from the edge) invades the 
region of interest. 

An immediate application of sCOLA is as a method for performing zoom-in simulations, 
which we illustrated in this paper. The fact that the N-body side of sCOLA is needed to 
perform calculations of the potential only for the near field, implies that sCOLA-based N- 
body codes can be made embarrassingly parallel, thus opening the door for efficiently tiling a 
volume of interest using grid computing in what may become Cosmology@Home. Moreover, 
combining sCOLA with the original COLA method can be useful for cheaply generating 
large ensembles of accurate mock halo catalogs required to study uncertainties of large-scale 
structure observables. Surveys with large aspect ratios, such as pencil-beam surveys, can 
especially benefit, as sCOLA can easily capture the effects of large-scale transverse modes 
without the need to substantially increase the simulated volume. 

pyCOLA can be found on the following URL: https: //bitbucket. org/tassev/pycola/. 
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A The details behind sCOLA 


In this appendix we write down the exact first equation of the sCOLA algorithm, given 
schematically by eq. (2.4). One should read this appendix as a continuation of Appendix 
A of TZE. The preliminaries, notation, and original temporal-domain COLA (tCOLA) are 
presented there. Here we present the modifications to the original tCOLA method needed to 
extend COLA to the spatial domain. 


A.l The sCOLA method 

The COLA method (both tCOLA and sCOLA) dictates that we integrate the equation of 
motion for CDM particles in a frame comoving with observers following trajectories specified 
by LPT. Thus, we solve for the residual displacement of CDM particles: 

s res = s - D1S1 - D 2 s 2 ■ (A.l) 

For tCOLA, the equation of motion is written as (see eq. A.6 in TZE) 

T 2 [s res ] = ~n M ad x d- 2 8(x,a ) - T 2 [D\}si - T 2 [D 2 ]s 2 . (A.2) 

For sCOLA, we rewrite it as 

T 2 [s ies ] = ~n M ad x dZ OLA ’- 2 5 (x, a) - T 2 ^]*? 01 ^ - T 2 [D 2 ]s% OLA . (A.3) 

Thus, the only difference between sCOLA and tCOLA is adding the COLA superscript 
where appropriate, reminding us that those quantities are calculated in the small COLA box, 
covering only a fraction of the original volume. The calculation of the LPT displacement 
fields, sfP LA , in the COLA box is presented in Appendix B. 

Below we present results for the combined tCOLA/sCOLA discretization as employed 
in pyCOLA. Therefore, as described in TZE, next we discretize the operator T only on the 
left-hand side of the above equation, using as a velocity variable v res = T[s res ], which is the 
CDM velocity as measured by an LPT observer. 


A.2 Discretizing T in sCOLA in the standard approach 

Similar to TZE, we discretize the sCOLA equations using a kick-drift-kick scheme. That 
requires writing down a kick operator K, a drift operator D, as well as introducing two more 
operators, L±, for converting between the standard particle velocity and the residual velocity 
as measured by LPT observers. 

In the “standard” prescription (see A.3.1 of TZE), the operators D and L± of tCOLA 
and sCOLA are identical. The kick operator in sCOLA, however, is modified (compare with 
eq. A.7 of TZE): 


D(a^,a/;a c ): x{ai) x(a f ) = 


E- {p>i 5 5 a c ) • 


a f 

(a/) = x(ai) + v(a c ) J - 


da 

Qifl) 


+ 


(A.4) 


+(Di(af) — Di(ai))s\ + ( D 2 (af ) — D 2 {ai))s 2 , 

' “/ 

v(ai) v(af) = v{ai) 


/ 


a/a c 


Q(a 


da 


x 


(— ln M a c d x d^ 01jA ’- 2 8(x,a c ) - - T 2 [D 2 }(a c )s^ OLA ^ , 
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where as in TZE we dropped the (res) subscript from v res for brevity. 

Time evolution between ao and a n + 1 in sCOLA is achieved in exactly the same way as 
in tCOLA, by applying the following operator on (®(ao), v(clq) = T[s](ao)): 

L+(a n+ i) ^Y\K(a i+1 / 2 , a i+1 - a i+ i)B(ai, a i+ r, a i+ 1 / 2 )K(ai, a i+1 / 2 ] L_(a 0 ) . (A.5) 

As in TZE, here we defined 

L±(a) : v(a) v(a) = v(a) ± ^T[Di](a)si + T[D2](a)s2^J , (A.6) 

which as mentioned above, first transforms the initial conditions for v to the rest frame of 
LPT observers (L_), and then adds back the LPT velocities at the end (L + ). 


A.3 Modified discretization of T in sCOLA 


As done in Appendix A.3.2 of TZE, here we present a modified discretization scheme of T, 
which allows one to improve the convergence properties of COLA. As above, the L± and D 
operators of tCOLA and sCOLA in this scheme are identical. The kick operator in sCOLA, 
however, is modified (compare with eq. A. 15 of TZE): 


D(a^,a/;a c ) : 


a f 

/x , x , x v ( a c ) f u ( a ) _ 
x ( ai ) x { CLf ) = x ( ai ) H--—- / —, x da + 


u(a c ) J Q(a ) 


(A.7) 


CLi 


+(Di(a,f) — Di(ai))si + ( D 2 (a,f ) — D 2 (ai))s 2 , 

u(a f ) - u(ai) 


K(a,,a/;a c ): v(a,i) i-» v(a f ) = v(a,i) + 


T[u](a c 


-^Ma c d x dn OLA ’ 2 5(x,a c ) -T 2 [Di](a c )sf OLA -T 2 [D 2 ](a c )s$ OLA \ , 


where we again dropped the (res) subscript from v res for brevity. The equations presented 
above give the evolution equations behind the sCOLA algorithm presented in this paper. 

However, in practice, we found that using the above prescription results in density fields 
from sCOLA which are artificially smoothed. The reason is due to a mismatch at very short 
scales between the LPT displacement fields calculated in the sCOLA box and the original 
simulation volume. Such a mismatch at deeply non-linear scales has no deep physical origin 
and we view it as an artifact of the different algorithms we use for calculating the LPT 
fields, combined with numerical precision issues. Removing that artifact is trivial: we boxcar 
average 9 (with width of 2 particle cells) all LPT displacement fields entering the operators 
L±, K and D. That eliminates most short-scale noise, which otherwise artificially smoothed 
our final results. We would like to emphasize that no information is lost in such an averaging 
as LPT carries little (if any) information on the scale of the averaging window we chose. 
Moreover, the convergence of our scheme to the “truth” when the number of timesteps and 
sCOLA volume are brought to infinity is not affected since in that case we would add and 
subtract the same averaged LPT displacements. Note that the averaging is done only on the 
displacements entering the evolution operators, and not on the initial conditions, as we want 
to keep the initial conditions intact. 

9 For different simulation set-ups, one may need to experiment with varying the boxcar window, or using a 
Gaussian averaging instead. 
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B Calculating LPT displacements using force evaluations 


In this appendix we write down the exact second equation of the sCOLA algorithm, given 
schematically by eq. (2.5), and describe how we solve it in pyCOLA. The non-schematic 
version of equation (2.5) reads: 

OLA + T 2 [D 2 ]s % ola = -^Q M aVVcoLA^[*LPT](*LPT) (solved in PT), 

(B.l) 


where x L pt = q + D\s i + D 2 s 2 . 

The usual way to derive the LPT displacement fields, s\ and s 2 , is using Fourier tech¬ 
niques [27]. However, the method we employ is based on [28], where the equation above is 
solved in real space by calculating the gravitational forces due to the density field entering 
above at early enough times that one can consider the results perturbative. The right-hand 
side above is completely determined by the initial conditions in the full volume. So, all we 
need to do is extract its first- and second-order pieces and match those to the corresponding 
terms on the left-hand side. We do that explicitly below, following [28]. This is advantageous 
in our case because the displacement fields in the equation above live in two separate boxes 
- the COLA box (s^ 2 >LA ) and the full box (sp 2 ). The results are summarized below, and 
are implemented in the public pyCOLA code. These results include modifications to [28] 
(allowing /3 = 1 in the notation below) to obtain high-precision LPT displacements for the 
sCOLA box. For further information, see [28]. 

pyCOLA offers a choice in the calculation of the first-order and second-order displace¬ 
ments, sp 2 >LA , in the COLA volume, as well as s 2 in the full volume at redshift zero: using 
a 2-pt or a 4-pt (denoted by superscript) rule. Depending on the rule, the first/second order 
displacements receive higher-order corrections as follows: at third/fourth order when using 
the 2pt rule; at fifth/sixth order when using the 4pt rule. The displacement fields are given 
by (77 below stands for the displacement field in the COLA box or the full box): 


where 


c ^,2pt 

s i 

c ^7,2pt 

S 2 


S 1 


S 2 


-^[F v {g,/3g 2 )-F v {-g,l3g 2 )\ , 

-^2 [F v {g,f3g 2 ) + F v (-g,(3g 2 )] , 

-4-2—T k(<?,/% 2 ) - F v (-g,Pg 2 ) - 
Zg a z — 1 

- 4 (-C; ( ag, /3a 2 g 2 ) - F v (- ag , /3a 2 g 2 ) ) , 

~^2^j[ F v( 9 ^9 2 ) + F v(-9,(3g 2 ) ~ 

~^ (-C; {ag, /3a 2 g 2 ) + F v {-ag, /3a 2 g 2 ) ) , 


a = (3/10)f ^ 143 if [3 = 1, 
a = (3/7)f ^ 143 if 0 = 0 , 


(B-2) 


(B.3) 
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and a / 1 and g are arbitrary constants. The equations above become exact in the limit 
of g — y 0, ag — y 0. When writing the equations above in the COLA box to second order 
(i.e. g = COLA), one must use /? = 1 in order to be consistent in PT. If one wants to use 
the above equations for calculating S 2 in the full box, then one must set /? = 0 since (3 = 1 
requires one to provide S2- The factors of O^ 143 (^m being the matter density today) are 
needed to rescale the second order displacements to matter domination and are correct to 
0(max(lO -4 , g 3 /143)) in ACDM. 

The force 1 ^( 51 , 52 ) corresponds to the right-hand side of eq. (B.l). It is evaluated 
numerically by using the force calculation function (using a 4-pt finite difference scheme) of 
pyCOLA on the PM grid. Analytically, F v is given by: 


Fr,{gi,g2) 


vv- 2 5 


Q + 5'i'Si + <?2 ^ m 1/143 s 2 


(B.4) 


where g\ and g 2 correspond to the first and second order growth factors, assumed to be deep 
in matter domination (see the arguments of F in eq. (B.2) above). 

The fractional overdensity, 5, in pyCOLA is calculated using cloud-in-cell assignment 
from a grid of particles (at Lagrangian positions q) displaced by the first and second order 
displacement fields: g\S\ + g2CtjJ^ US $ 2 - Note that implicitly for each particle at Lagrangian 
position g, the force F rj (gi, ^ 2 ) is evaluated at its corresponding Eulerian position: q + g\S\ + 
92^m^ U3s 2- The choice of 77 (77 = COLA or 77 = “full”) determines which particles (only 
those in the COLA volume, or all particles in the full volume) are included in the calculation 
of the potential, V“ 2 5. 
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